#
# pm-svr.py
# author:           Viet Hung
# description:      PM regression using Support Vector Regression/MLR
#

import os
import sys
import numpy
import math
from sklearn import svm, linear_model, ensemble, tree
from sklearn.metrics import mean_squared_error
from math import sqrt
import pickle
import datetime
from datetime import timedelta
import pylab as pl
from scipy import stats

#
# Get the feature base on the feature_mask
#
# Input:
#   data            all data
#   feature_mask    the feature mask
#
# Output:
#   feature_vectors     a list of feature vector using the feature mask
#                       
def get_feature_vector(data, mask):
    row = len(data)
    col = 0

    new_data = []

    for row in data:
        a = []
        for i in range(len(mask)):
            if mask[i] == 1:
                a += [row[i]]

        new_data += [a]

    return new_data

#
# Use svr regression for predicting pm value
#
# Input:
#   training_data       training data
#   test_data           test data
#   mask                the feature mask
#   index               the index of PM value (3, 4, 5)
#
# Output:
#   r_square            the square of correlation error
#   rms                 the root means square errors
#
def svr_regression(training_data, test_data, mask, index = 0):
    # Get the training, test parameters
    n_training = len(training_data)
    training_feature = get_feature_vector(training_data, mask)
    training_target = training_data[:, index].transpose().tolist()

    n_testing = len(test_data)
    test_feature = get_feature_vector(test_data, mask)
    test_target = test_data[:, index].transpose().tolist()

    # training_feature, training_target, test_feature, test_target = normalized(training_feature, training_target, test_feature, test_target)
     
    # Train svr model
    f = linear_model.LinearRegression()
    # f = svm.SVR(kernel="poly", degree=3, epsilon=2.0, C = 2.0, gamma = 0.75, coef0 = 2.5)
    f.fit(training_feature, training_target)
    
    print f
    print "intercept", f.intercept_, "coef", f.coef_

    r_square, rms, re = calc_error(f, test_feature, test_target)
    
    return r_square, rms, re

#
# Compute the predict error
#
# Input:
#       f           regression model
#       X           feature vectors
#       y           targets
#
def calc_error(f, X, y):
    if len(X) != len(y): 
        return -1, -1, -1

    y1 = numpy.array(f.predict(X))

    # 
    # slope, intercept, r, p, std_err = stats.linregress(y, y1)
    # print "R-Squared", r**2, "STD-Err", std_err

    # calculate error
    r = numpy.corrcoef([y, y1])
    rms = sqrt(mean_squared_error(y, y1))
    re = 0
    for i in range(len(y)):
        re += abs((y[i] - y1[i]) / y[i])

    return r[0, 1]*r[0, 1], rms, re / len(y) * 100.0

def transform(x, a, b):
    return [a*x[i] + b for i in range(len(x))]

#
# Data normalization: scale all the features and targets into the range [-1 1]
#
def normalized(training_feature, training_target, test_feature, test_target):
    n_feature = len(training_feature[1])

    max_val = [0. for i in range(n_feature)]
    min_val = [0. for i in range(n_feature)]
    mid_val = [0. for i in range(n_feature)]
    val_range = [0. for i in range(n_feature)]
	
    tmp1 = numpy.array(training_feature)
    tmp2 = numpy.array(test_feature)
    for i in range(1, n_feature):
        tmp3 = numpy.hstack((tmp1[:, i], tmp2[:, i]))
        max_val[i] = tmp3.max()
        min_val[i] = tmp3.min()
        val_range[i] = max_val[i] - min_val[i]
        mid_val[i] = (max_val[i] + min_val[i])/2

    #print training_target
    tmp1 = numpy.hstack((numpy.array(training_target),numpy.array(test_target)))
    max_target = tmp1.max()
    min_target = tmp1.min()
    target_range = max_target - min_target
    mid_target = (max_target + min_target)/2

    for i in range(len(training_feature)):
        for j in range(1, n_feature):
            training_feature[i][j] = (training_feature[i][j] - mid_val[j])/val_range[j] 

        #training_target[i] = math.log(training_target[i])

    for i in range(len(test_feature)):
        for j in range(1, n_feature):
            test_feature[i][j] = (test_feature[i][j] - mid_val[j])/val_range[j]
        
        #test_target[i] = math.log(test_target[i])
   
    return training_feature, training_target, test_feature, test_target

def add_years(d, years):
    """Return a date that's 'years' years after the date (or datetime)
    object 'd'. Return the same calendar date (month and day) in the
    destination year, if it exists, otherwise use the following day
    (thus changing February 29 to March 1).
    """
    try:
        return d.replace(year = d.year + years)
    except ValueError:
        return d + (date(d.year + years, 1, 1) - date(d.year, 1, 1))

def regression_by_month(FILE):
    # Use for selecting feature for training
    FEATURE_MASK = [0, 0, 0, 0, 1, 1]
    data = numpy.loadtxt(FILE, delimiter = ",", usecols=(1,5,6,7,8,9))
    
    begin_time = 2012.0
    end_time = 2012.5
    for i in range(5):
        print "MLR using testing data in year", begin_time

        testing_data = numpy.array(filter(lambda x: x[0] >= begin_time and x[0] < end_time, data))
        training_data = numpy.array(filter(lambda x: x[0] < begin_time or x[0] >= end_time, data))

        print "Training data", len(training_data), "rows"
        print "Testing data", len(testing_data), "rows"
        if len(training_data) > 0 and len(testing_data) > 0:
            r_square, rms, re = svr_regression(training_data, testing_data, FEATURE_MASK, index = 2)
            print "R2, RMS, RE = (%f, %f, %f) " % (r_square, rms, re)
        print "\n"
        print "-------------------------------------------------------------------------------"

        begin_time += 0.5
        end_time += 0.5

def regression_by_year(FILE):
    # Use for selecting feature for training
    FEATURE_MASK = [0, 0, 0, 0, 1, 1]
    data = numpy.loadtxt(FILE, delimiter = ",", usecols=(1,5,6,7,8,9))
    
    begin_time = 2011
    end_time = 2012
    for i in range(3):
        begin_time += 1
        end_time += 1
        print "MLR using testing data in year", begin_time

        testing_data = numpy.array(filter(lambda x: x[0] >= begin_time and x[0] < end_time, data))
        training_data = numpy.array(filter(lambda x: x[0] < begin_time or x[0] >= end_time, data))

        print "Training data", len(training_data), "rows"
        print "Testing data", len(testing_data), "rows"
        if len(training_data) > 0 and len(testing_data) > 0:
            r_square, rms, re = svr_regression(training_data, testing_data, FEATURE_MASK, index = 2)
            print "R2, RMS, RE = (%f, %f, %f) " % (r_square, rms, re)
        print "\n"
        print "-------------------------------------------------------------------------------"
def regression_by_year_2(FILE):
    # Use for selecting feature for training
    FEATURE_MASK = [0, 0, 0, 0, 1, 1, 1]
    data = numpy.loadtxt(FILE, delimiter = ",", usecols=(1,5,6,7,8,9,11))
    
    begin_time = 2011
    end_time = 2012
    for i in range(3):
        begin_time += 1
        end_time += 1
        print "MLR using testing data in year", begin_time

        testing_data = numpy.array(filter(lambda x: x[0] >= begin_time and x[0] < end_time, data))
        training_data = numpy.array(filter(lambda x: x[0] < begin_time or x[0] >= end_time, data))

        print "Training data", len(training_data), "rows"
        print "Testing data", len(testing_data), "rows"
        if len(training_data) > 0 and len(testing_data) > 0:
            r_square, rms, re = svr_regression(training_data, testing_data, FEATURE_MASK, index = 2)
            print "R2, RMS, RE = (%f, %f, %f) " % (r_square, rms, re)
        print "\n"
        print "-------------------------------------------------------------------------------"

def regression_by_station(FILE):
    # Use for selecting feature for training
    FEATURE_MASK = [0, 0, 0, 0, 1, 1]
    data = numpy.loadtxt(FILE, delimiter = ",", usecols=(0,5,6,7,8,9))
    
    for i in range(1, 6):
        print "MLR using testing data in station", i

        testing_data = numpy.array(filter(lambda x: x[0] == i, data))
        training_data = numpy.array(filter(lambda x: x[0] != i, data))

        print "Training data", len(training_data), "rows"
        print "Testing data", len(testing_data), "rows"
        if len(training_data) > 0 and len(testing_data) > 0:
            r_square, rms, re = svr_regression(training_data, testing_data, FEATURE_MASK, index = 2)
            print "R2, RMS, RE = (%f, %f, %f) " % (r_square, rms, re)
        print "\n"
        print "-------------------------------------------------------------------------------"
def regression_by_station_2(FILE):
    # Use for selecting feature for training
    FEATURE_MASK = [0, 0, 0, 0, 1, 1,1]
    data = numpy.loadtxt(FILE, delimiter = ",", usecols=(0,5,6,7,8,9,11))
    
    for i in range(1, 6):
        print "MLR using testing data in station", i

        testing_data = numpy.array(filter(lambda x: x[0] == i, data))
        training_data = numpy.array(filter(lambda x: x[0] != i, data))

        print "Training data", len(training_data), "rows"
        print "Testing data", len(testing_data), "rows"
        if len(training_data) > 0 and len(testing_data) > 0:
            r_square, rms, re = svr_regression(training_data, testing_data, FEATURE_MASK, index = 2)
            print "R2, RMS, RE = (%f, %f, %f) " % (r_square, rms, re)
        print "\n"
        print "-------------------------------------------------------------------------------"
def regression_all_data(FILE):
    # Use for selecting feature for training
    FEATURE_MASK = [0, 0, 0, 0, 1, 1]
    data = numpy.loadtxt(FILE, delimiter = ",", usecols=(0,5,6,7,8,9))
    print "Training data", len(data), "rows"
    if len(data) > 0:
		# index: vi tri cua cot du lieu muon suy ra
        r_square, rms, re = svr_regression(data, data, FEATURE_MASK, index = 2)
        print "R2, RMS, RE = (%f, %f, %f) " % (r_square, rms, re)
    print "\n"
    print "-------------------------------------------------------------------------------"
def regression_all_data_2(FILE):
    # Use for selecting feature for training
    FEATURE_MASK = [0, 0, 0, 0, 1, 1,1]
    data = numpy.loadtxt(FILE, delimiter = ",", usecols=(0,5,6,7,8,9,11))
    print "Training data", len(data), "rows"
    if len(data) > 0:
        # index: vi tri cua cot du lieu muon suy ra
        r_square, rms, re = svr_regression(data, data, FEATURE_MASK, index = 2)
        print "R2, RMS, RE = (%f, %f, %f) " % (r_square, rms, re)
    print "\n"
    print "-------------------------------------------------------------------------------"
if __name__=="__main__":


    # MOD04 all
    #regression_all_data_2("data/mod.csv")
    #regression_by_station_2("data/mod.csv")
    #regression_by_year_2("data/mod.csv")
    
   